Vortices in Ginzburg-Landau billiards 



E. Akkermans and K. MallickQ 
Department of Physics, Technion, 32000 Haifa, Israel 

We present an analysis of the Ginzburg-Landau equations for the description of a two-dimensional 
superconductor in a bounded domain. Using the properties of a special integrability point of these 
equations which allows vortex solutions, we obtain a closed expression for the energy of the su- 
perconductor. The role of the boundary of the system is to provide a selection mechanism for the 
number of vortices. A geometrical interpretation of these results is presented and they are applied 
to the analysis of the magnetization recently measured on small superconducting disks. Problems 
related to the interaction and nucleation of vortices are discussed. 
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I. INTRODUCTION 



This article is devoted to the study of the existence and stability of vortex solutions of the Ginzburg-Landau 
equations in finite two-dimensional domains with boundaries. For the case of the infinite plane, a large amount of 
work has been devoted to such questions, motivated by the ubiquitous character of the Ginzburg-Landau equations 
either in the study of superconductors or for the Abelian Higgs model. It is known that for the infinite system 
the Ginzburg-Landau energy functional has a lower bound which is saturated for a special choice of the parameters 
(which corresponds to the limiting case between type I and type II superconductors). The minimizing equations 
at this special point, also called the dual point, are first order differential equations and they admit solutions with 
vortices. These different vortex solutions are classified by an integer of topological origin, namely the winding number 
of the order parameter which, roughly speaking, counts the number of vortices. 

Although the Ginzburg-Landau equations in finite domains have been investigated there is no generalization 
so far of the results concerning the dual point for finite two-dimensional domains (that we shall also call hereafter 
'billiards'). In this paper, we study such a generalization; besides a theoretical interest, our motivation has been 
triggered by a set of new experimental results obtained on small aluminium disks j^] in a mesoscopic regime where 
the radius R is comparable with both the coherence length £ and the London penetration length A. There, the 
magnetization, as a function of the applied magnetic field, presents a series of jumps with an overall shape reminiscent 
of type-II superconductors. Such a behaviour is very different from that of macroscopic systems: in the large, 
aluminium is a genuine type-I superconductor. These experiments were first analyzed in the framework of the 
linearized Ginzburg-Landau equations |J. Although this approach explains some of the observed features, it fails to 
provide a satisfactory quantitative agreement, and can not account for the large spatial variations of the magnetic 
field inside the sample (e.g. vortices). Numerical solutions |7) of the Ginzburg-Landau equations give a much better 
description of this phenomenon, and emphasize the importance of the nonlinear term. However, in this context, no 
analytical results are available for the full Ginzburg-Landau equations. Using our results on the dual point, we shall 
derive an analytical expression for the free energy and for the magnetization of a mesoscopic superconductor as a 
function of the applied magnetic field. 

The plan of this paper is as follows. In the section 2, we shall describe the main features and the known results 
about the dual point and the lower energy bound for an infinite system. Then, in section 3, we shall study the case of 
a finite domain and derive an expression of the free energy which, in addition to the previous lower bound, includes 
a contribution of the boundary. Section 4 contains a geometrical interpretation of these results. In section 5, our 
results are applied to analyze experimental data on small superconducting disks. In the conclusion, we propose some 
extensions of our work and a scenario for the nucleation of vortices at the boundary. 



II. THE DUAL POINT OF GINZBURG-LANDAU EQUATIONS IN AN INFINITE SYSTEM 

From now on, we study a superconducting billiard within the framework of the Ginzburg-Landau equations (this 
assumes that both the order parameter and the vector potential have a slow spatial variation) . The expression of the 
Ginzburg-Landau energy density a is given by 

2e - B 2 

a = a Q + a 2 \iP\ 2 + a 4 \^\ 4 + Oi|(V - i--A)i>\ 2 + — (1) 
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where tp = \ip\e lx is the complex-valued order parameter, B is the magnetic field and the aj's are real parameters. 
The coherence len 

A 2 = -t^(^) 2 ? 4 i , so that we obtain for the dimensionless free energy T 



The coherence length and the London penetration length are related to these parameters as follows || : £ 2 = t^t and 



T = I l\B\ 2 + K 2 \1- HT + |(V - iA)^\ 2 (2) 




where ip is measured in units of ipo = y \£2± ; B in units of j^p- , where 4>q = ^| is the quantum of flux. The lengths 
are measured in units of (the numerical factor V2 is for further convenience). The ratio k = ^ is the only 
free parameter in (2) and it determines, in the limit of an infinite system, whether the sample is a type-I or type-II 
superconductor ||. The integral is over the two-dimensional domain of the superconducting sample. 

The Ginzburg-Landau equations for the order parameter tp and for the magnetic field B = V x A are obtained from 
a variation of T: 
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where the current density j is given by 



(V - iA) 2 t/j = 2k 2 iI>(1 - \4>\ 2 ) 

V x B = 2j (3) 

j = Im(tp*Vtp) - \ip\ 2 A (4) 



The Ginzburg-Landau equations are nonlinear second order differential equations and their solutions are not known 
except for some particular cases. However, for the special value k = -j=, the equations for tp and A can be reduced to 
first order differential equations. This special point was first used by Sarma in his discussion of type-I vs type-II 
superconductors and then identified by Bogomol'nyi |IJ in the more general context of stability and integrability of 
classical solutions of some quantum field theories. We now review some of these properties of the Ginzburg-Landau 
free energy at the dual point. We shall use the following identity true for two dimensional systems 

|(V-L4)V| 2 = \T>iP\ 2 +V x j + B\tp\ 2 (5) 

where J*is the current density and the operator D is defined as T> = d x +id y — i(A x + iA y ). At the dual point, n = 
the expression (2) for T can be rewritten using the identity (||) as follows 

T= (\\B-l+ H 2 | 2 + \V^\ 2 + / (j + A).dl (6) 
Jn z Jan 

where the last integral over the boundary d£l of the system results from Stokes theorem. 



For an infinite system we impose, as in 
j — > at infinity. The boundary term in (j 



, that the system is superconducting at large distance, i.e. \tp\ — ► 1 and 
) then becomes 



(j + A).dl=<p (-±-+A).dl (7) 



_ J 
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This last integral is known as the London fluxoid. It is quantized and using (^) one shows that it is equal to 
§an VX'dH = 27m, where \ is the phase of the order parameter. The integer n is the winding number of the order 
parameter tp and as such is a topological characteristic of the system. The extremal values of J- , namely T = 27m, are 
obtained when the bulk integral in (^) vanishes identically, i.e. when the two Bogomol'nyi Q equations are satisfied: 

x>V = o 

B = 1-\tP\ 2 (8) 
These two equations can be decoupled and one obtains that \ip\ is a solution of the second order nonlinear equation 

V 2 HVf = 2(|V| 2 -l) (9) 

This equation is related to the Liouville equation. The set of equations (|^J^) has been obtained without any assumption 
on the nature of the magnetic field and appears in various other situations, e.g. Higgs ||, Yang-Mills || and Chern- 
Simons pj field theories. It was proven that these equations admit families of vortex solutions For infinite 
systems, it can be shown that each vortex carries one flux quantum and that the winding number n is equal to the 
number of vortices in the system. However, for an infinite system there is no mechanism to select the value of n, 
which only plays the role of a classifying parameter. It will be precisely the role of the boundary of a finite system to 
introduce such a selection mechanism and to determine n, according to the applied magnetic field. 



III. THE FINITE SIZE SYSTEM 



From now on, we shall study finite size systems in an external magnetic field. The question then arises to know if 
such systems can sustain stable vortex solutions and how they behave as a function of the applied field. A simplified 
version, without applied magnetic field, was extensively studied by Bethuel et al. Q). In this work the mechanism 
for vortex creation is based on Dirichlct boundary conditions of the type ip = f on <9f2 and where / is a complex 
function of degree n. In the London limit, namely k — > oo, \ip\ is 1 almost everywhere but, because of the degree 
n imposed on the boundary, ip must vanish n times in the bulk therefore leading to vortices. Moreover, numerical 
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simulations of the Ginzburg-Landau equations for a long and thin parallepiped in a uniform magnetic field Jil| show 
the existence of stationary vortex solutions whose number depends on the applied magnetic field. These simulations 
then indicate that the physical picture derived for k = remains valid for quite a large range of values of k, and 

that the corresponding change of free energy is small (see also E^] ) . Indeed, for finite systems their size R becomes 
relevant as a new length scale so that k is no more the only dimensionless parameter, and its value does not alone 
control the physics at the mesoscopic scale as it does for infinite systems. We shall study the case n — -?=, i.e. the 
dual point, and extend the previous approach to a system with finite size where boundary effects are important. 

In a finite system, there are in general non-zero edge currents and the order parameter is not equal to 1 on the 
boundary. Hence, the identification of the boundary integral in (||) with the fluxoid (0) is not possible anymore, and 
the free energy can not be minimized just by imposing Bogomol'nyi equations (||). However, the currents on the 
boundary of the system screen the external magnetic field and therefore produce a magnetic moment (a circulation) 
opposite to the direction of the field, whereas vortices in the bulk of the system produce a magnetic moment along the 
direction of the applied field. Hence currents in the bulk circulate in a direction opposite to those at the boundary. 
If one assumes cylindrical symmetry, the current density Jhas only an azimuthal component, with opposite signs 
in the bulk and on the edge of the system (the radial component of J* is zero since J is divergenceless). Thus, there 
exists a circle T on which / vanishes [ fl3|| . This allows us to separate the domain fl into two concentric subdomains 
fi = Oi U O2 such that the boundary dili is the curve T (see Fig. I). On d£l\, the current density Jis zero, therefore 
one has 

f.dl=(f rj^.dl = (10) 

Thus one deduces as above that Bogomol'nyi and Liouvillc equations are valid in the finite domain f^i as in the case 
of the infinite plane. The existence of vortices in a finite domain such as Oi was checked using a numerical solution 
JLlj of (^): assuming cylindrical symmetry, one shows that \ip\ vanishes as a power law at the center of the disk, hence 
there is a vortex in the center (more precisely a multi- vortex whose multiplicity is determined by the exponent of the 
power law); moreover, saturates very rapidly to a constant value close to one for lengths larger than A. The same 
conclusion can be reached by defining f(r) — —ln\ip\ 2 and linearizing (^) around \tp\ — 1. Then, / satisfies V 2 / = 2/ 
whose general solution is f(r) = aIo(r^/~2) + 6_ftTo(rv / 2). From the behaviour of the Bessel functions Iq and Kq, one 
obtains that for small r, vanishes as a power law and saturates rapidly to a constant of order one in a finite range 
of r's (in units of \V2). 

The magnetic flux <J?(f2i) through is calculated, in units of the flux quantum </> , using the fluxoid and ( |l0| ) so 
that 



*(n x ) = n-/ £f 



As before n is the winding number, i.e. §qq Vx.dZ = 27m, as well as the number of vortices. The free energy in fli is 

F(Oi) = 27m. (II) 

We consider now the contribution of fl 2 to the free energy. It is given by (2) and can be rewritten using the phase 
and the modulus of the order parameter ip, as 

T(n 2 ) = f (vh) 2 + H 2 |v x -A 2 + ^+ (1 " if |2)2 (12) 

We know, from the London equation, that both the magnetic field and the vector potential decrease rapidly (exponen- 
tially) away from the boundary Oil of the system over a distance of order 1 in units of X\^2. Over the same distance, 
at the dual point, saturates to unity. One can thus estimate the integral (|l^) using an elementary version of the 
saddle-point method. We assume cylindrical symmetry, and we neglect the term (V|i/>|) 2 on the boundary because 
of the covariant Neumann boundary conditions at the interface between a superconductor and an insulator ||. We 
obtain for the free energy the expression 

F{<h)*f H 2 |V X -Al 2 + ^+ (1 -y |2)2 (13) 



no 
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where the integral is now over the boundary of the system. To go further, we need to implement boundary conditions 
for the magnetic field B(R) and the vector potential A(R). The choice B(R) = B e , where B e is the external imposed 
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field, corresponds to the geometry of an infinitely long cylinder, where the flux lines are not distorted outside the 
system. This boundary condition is not adapted to describe a flat thin disk. A more suitable choice is provided by 
demanding <fi = </> e , which means that the total flux through the disk is identical to the flux of the external field, 
although flux lines are distorted by the superconducting sample. One can check, using numerical simulations 
that this condition is well satisfied if the disk is thin enough. The boundary condition <f> = 4> e implies that the vector 
potential is identified by continuity to its external applied value A e which has only the azimuthal component 
It should be noticed that the magnetic field B has a non monotonous variation: it is low in the bulk, larger than 
B e near the edge of the system, because of the distortion of flux lines, and eventually equal to its applied value far 
outside the system 

Different choices of boundary conditions will give rise to different limits for a very large system (i.e. for R — ► oo). 
The limit of an infinitely long cylinder corresponds to a superconducting bulk sample, whereas the limit of a thin 
disk, which is the case we consider here, corresponds to a superconducting thin sheet. 

The formula (|l3|) for the free energy is similar to the Little-Parks expression || for a quasi one-dimensional hollow 
cylinder in a uniform applied magnetic field. The minimization with respect to \ip\ gives 1 — \ip\ 2 — |Vx ~ A\ 2 , such 



that (13) can be written as 



nv.,)=f \V X -A e \ 2 + \B 2 -±\V X -A e \ 4 (14) 



IdQ. 

Performing the integral over the boundary of the system, we obtain 

1 rtn s , n2 I.Aa/2 3 



27r~ m) = ~ir {n - 0e) " 2 ( ^r } (n _ 0e) (15) 



We have neglected the contribution of the B 2 term, which is similar to the first term in the r.h.s. of ([15]) but smaller 
by a factor of the order (X/R) 2 . The integer n which appears in (15) is the same as in (fTj), since the order parameter 



ij) is the same function in both subdomains. The circulation of its phase x (the winding number) counts the number of 
zeroes in the domain f2, i.e. the number of vortices. The thermodynamic Gibbs potential Q of the system is obtained 
from J- (fti) + T{VL2) by a Legendre transformation so that 

±g(n, =n+^(n- e ) 2 - \{^?{n - e ) 4 - 2 ^ 2 (16) 

This relation consists in a set of quartic functions indexed by the integer n. The minimum of the Gibbs potential is the 
envelop curve defined by the equation ^\(p c = 0, i.e. the system chooses its winding number n in order to minimize 
Q. This provides a relation between the number n of vortices in the system and the applied magnetic field <f> e . In the 
limit of a large enough y, the quartic term is negligible and the Gibbs potential reduces to a set of parabolas (Fig. 2). 
The winding number n is then given by the integer part 

The magnetization M — —■§§-> of the system, is given by 

-M = ^(^_n)-^ e (18) 

For (f> e smaller that , we have n = and (— M) increases linearly with the external flux. This corresponds to 
the London regime before the first vortex enters the system. The field Hi at which the first vortex enters the system 
corresponds to Q(n = 0) = Q(n = 1), i.e. to 

Hl -^VM~x + ^ (19) 

The subsequent vortices enter one by one for each crossing Q(n+ 1) = Q(n); this happens periodically in the applied 
field, with a period equal to 

AH=^ (20) 



5 



This gives rise to a discontinuity of the magnetization AM = 2 ^ A . 

There is a qualitative similarity between the results we derived using the Bogomol'nyi equations within the domain 
Sli and those obtained from a linearised version of the Ginzburg-Landau functional ||. But the two approaches differ 
in their quantitative predictions due to the importance of the non linear term. An illustration of this is given in the 
next section. 



IV. A GEOMETRICAL INTERPRETATION 



In an infinite system, one shows using the boundary conditions \ip\ — ► 1 and j — ► at infinity that equation (Q) 
implies the quantization of the magnetic flux: 



B.dS = n (21) 

R 2 

If one interprets B as a curvature, this relation is analogous to the Gauss-Bonnet theorem, which states that, for a 
compact manifold M without boundary, the integral of the Gaussian curvature K over the surface is equal to the 
Euler-Poincare characterics x °f the manifold, which is a topological invariant integer: 

K = X (22) 

M 

This result is in fact more than an analogy and at the dual point the Ginzburg-Landau functional has a useful 



geometrical interpretation that we shall now highlight |17 



The Ginzburg-Landau functional (2) for the energy of a superconductor corresponds to a U(l) gauge symmetry. 
For that symmetry, one can identify an abelian one-form connexion and a two-form curvature f2 respectively given 
by the vector potential A and by the magnetic field B. The complex order parameter ip is a section of the U(l) 
fiber and the basis manifold is the infinite plane. The above boundary conditions allow to map the plan onto a 
compact manifold without boundary, namely the sphere. Topological invariants of the problem are obtained from the 
Chern-Weil invariant polynomial P(f2) = det(l + -^Q) (see e.g. |Q) For the U(l) bundle over a sphere there is only 
one Chern class, f2(= B). The integral of that Chern class over the basis manifold is a topological invariant integer 
called a Chern number and is precisely given by (^lj). This Chern number plays a role similar to %■ At the dual point, 
the energy is T = J B = 27m; this fact can be translated in geometrical terms by saying that the extremal free energy 
is identical to a topological invariant of the problem namely its Chern number. 

When the basis manifold M has a boundary dM which is not a geodesic, the integral of its curvature is neither 
an integer nor a topological invariant. The Gauss-Bonnet theorem is then generalized by adding a boundary term so 
that the Euler-Poincare characteristics x is now given by the relation 

!/■„,„ 1 



X = — / KdS+—i k g dl (23) 
2-kJm 2tt J aM 

where K and k g are respectively the Gaussian curvature of the manifold and the geodesic curvature of the boundary. 
A similar result holds for the U(l) problem in a bounded domain (or in an infinite domain with boundary conditions 
different from those chosen above). In that case the relation equivalent to ( |23| ) is given by the fluxoid relation 



As before, B is the curvature of the connexion and the current density plays here the role of a geodesic curvature 
p7[ . The expression obtained in (|l^) for the Bogomol'nyi free energy of a system with a boundary can be rewritten 



(25) 



The boundary correction is a function 7/ of the geodesic curvature. The results obtained in the preceeding section 
show that the geodesic curvature is given by n — 4> e for a cylindrically symmetric system and using the appropriate 
approximations we determined the function r\ as being an even fourth order polynomial in the geodesic curvature. 

This geometric interpretation makes us believe that an expression such as (^5|) is fairly general. It could be well 
suited, as an Ansatz, to describe finite systems which are known to have a topological description in the infinite limit 



(for instance, a suitable generalisation of (25) to other symmetry group like SU(2) could describe some phases of 
superfluids in a bounded domain). 
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V. COMPARISON WITH THE EXPERIMENTAL RESULTS 



In order to compare our results to the experimental data given in |7| , we shall consider the limit where the radius R 
is larger than A and £, (typically, R ~ 10A is considered experimentally (HQ]). The thickness of the system is assumed 
to be small enough so that we can neglect variations of both the magnetic field and the order parameter along 
the cylinder axis. We shall prove that within these approximations, the expression ( |l6| ) captures the main features 
observed experimentally i.e. the behaviour of the magnetization at low fields (before the first discontinuity), the 
periodicity and the linear behaviour between the successive jumps and provides a fairly good quantitative agreement. 

The magnetization curve in fig. 3 (plotted for the ratio = 0.14, chosen arbitrarily) agrees qualitatively with both 
the numerical and experimental curves of fig. 3 in M. Besides, taking the experimental parameters of [Q, namely 
R = 1.2/xm and A(T) = 84nm at T = OAK, we compute from our expressions Hi = 25G and AH = 4.6G. These 
values agree with the results of [Q to within a few percent. We emphasize that H\ scales like j=r, whereas AH scales 
like -^j- in accordance with We calculate the ratio of the magnetization jumps to the maximum value of M to be 
0.20 as compared to a measured value of 0.22. The total number of jumps scales like R 2 and the upper critical field 
is independent of R in our theory. This to agrees fairly well with the experimental data ||[7| . We emphasize that this 
quantitative agreement cannot be obtained from the linearized Ginzburg-Landau equations. 



VI. CONCLUSION 



We have investigated the problems of the existence and the stability of vortices in a two-dimensional bounded 
superconducting system. To that purpose, we have used the Ginzburg-Landau energy functional at the special dual 
point characterized by the value k = which corresponds, for an infinite system, to a superconductor between type 
I and type II. We have shown for that case that it is still possible to obtain vortex solutions at the thermodynamic 
equilibrium. In contrast to the Bogomol'nyi solution obtained for the infinite plane, where there is no definite value 
for the number n of vortices, there is a selection mechanism for a finite billiard in an external magnetic field that 
allows to compute the number of vortices as a function of the applied field. Our reasoning is based on the existence 
of a contour T which allows to separate the system in two parts, the bulk containing the vortices and the edge with 
screening currents. As such the system might be viewed as a kind of two-dimensional Josephson junction or weak 
superfluid link flS|| . Vortices enter (or are expelled from) the superconductor through T, if the applied magnetic field 
is increased (or decreased). 

Although we considered the special case of the dual point, our analysis provides a satisfactory quantitative descrip- 
tion of the behaviour observed experimentally on such small superconducting billiards at least for the regime where 
the density of vortices is low enough i.e. for small applied magnetic fields. For higher fields, the expression ( |l6| ) does 
not describe properly the tail of the magnetization curve where both the periodicity AH and the amplitude AM of 
the jumps decrease and eventually vanish. In this regime, the vortices interact both between themselves and with the 
edge currents. This regime can be studied by a perturbative analysis around the dual point [ pp| . 

So far, we have studied only equilibrium states. On the basis of numerical simulations of the time-dependent 
Ginzburg-Landau equations 1 1 1 , we propose a mechanism for vortex nucleation. At each discontinuity of the winding 
number, namely when (f> e — is a half integer, there is a nodal line, joining the center of the system to its boundary, 
along which the order parameter ip vanishes and where its phase is ill defined. This might be interpreted as an opening 
of the ring of the screening currents which allows a flux line of the external field to enter the system. In this case, 
we expect the contour T to coincide with the boundary of the system. The existence of such a nodal line has been 
discussed in the context of Ginzburg-Landau equations and for the related Aharonov-Bohm problem of a magnetic 
flux line piercing either the infinite plane pl[ or a finite domain |^,^2| . 
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FIG. 1. Schematic set-up of the total system with the two subdomains £7i and f^2 separated by the contour V. 
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Applied magnetic flux 



FIG. 2. T{n,4»e) given by (1) plotted as a function of the applied magnetic flux <f) e for various values of the integer n and 
for X/R — 0.14. The free energy is the envelop of the ensemble of parabolas. 
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Applied magnetic flux 



FIG. 3. The magnetization — M as given by (15) as a function of the applied magnetic flux 4> e and for X/R = 0.14 . 
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